Take-home Exercise 2

Spatio-temporal Analysis of COVID-19 Vaccination Trends at the Sub-district Level, DKI Jakarta

Published

February 26, 2023

Modified

March 1, 2023

1 Context

1.1 Problem Statement

Understanding how Covid-19 vaccination rate is distributed around DKI Jakarta and how it changes over time.

1.2 Data

Data Format (Type) Description Source
Riwayat File Vaksinasi DKI Jakarta Excel (Aspatial) Details the daily number of vaccinations in Jakarta by person type and sub-district (Only the first day of each month is used) Riwatat File Vaksinasi DKI Jakarta (2022)
DKI Jakarta Adminstration Boundary 2019 Shapefile (Geospatial) DKI Jakarta administrative boundaries Indonesia Geospasial (2019)
Note

For the Vaksinasi DKI Jakarta data for March 2022, there was no available link to 1 March 2022 data. Thus for March 2022, the first day available, 2 March 2022, will be used instead.

2 Setting Up

2.1 Installing and Loading Packages

pacman::p_load(sf, tmap, sfdep, tidyverse, zoo, readxl, plotly, Kendall, ggplot2, dplyr)

  • sf: importing and handling geospatial data

  • tidyverse: wrangling attribute data

  • sfdep: newer package of spdep (does not require conversion to sp)

  • tmap: prepare cartographic quality chropleth maps

  • zoo: for regular and irregular time series

  • readxl: read excel (xlsx) format

  • plotly: interactive web graphics from ggplot2 graphs

  • Kendall: to perform Mann-Kendall test

  • ggplot2: plot temporal trend

  • dplyr: join_by() function

2.2 DKI Jakarta Administrative Boundary 2019

2.2.1 Importing Data

To import a shapefile, use st_read() where dsn is the targeted directory and layer is the name of the files.

jkt <- st_read(
  dsn = "data/geospatial",
  layer = "BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA"
)
Reading layer `BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA' from data source 
  `D:\bellekwang\IS415\take_home_ex\THE2\data\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 269 features and 161 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 106.3831 ymin: -6.370815 xmax: 106.9728 ymax: -5.184322
Geodetic CRS:  WGS 84
  • jkt is a multipolygon sf dataframe with 269 features

2.2.2 Projection

The current geodetic CRS is WGS 84 which is wrong. Using st_transform(), change the CRS from WGS 84 to DGN95 (EPSG:23845), the national projected coordinates systems of Indonesia TM-3 zone 54.1.

jkt <- st_transform(jkt, crs = 23845)
st_crs(jkt)
Coordinate Reference System:
  User input: EPSG:23845 
  wkt:
PROJCRS["DGN95 / Indonesia TM-3 zone 54.1",
    BASEGEOGCRS["DGN95",
        DATUM["Datum Geodesi Nasional 1995",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4755]],
    CONVERSION["Indonesia TM-3 zone 54.1",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",139.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9999,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",200000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",1500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre."],
        AREA["Indonesia - onshore east of 138°E."],
        BBOX[-9.19,138,-1.49,141.01]],
    ID["EPSG",23845]]

2.2.3 Data Cleaning

2.2.3.1 Outer Island

First, using tmap functions tm_shape() and tm_polygons(), plot out jkt to understand the area covered. I will also be adding the “DESA” label to identify which are the areas to exclude.

tmap_mode("view")
tm_shape(jkt) +
  tm_polygons() +
  tm_text("DESA") +
tm_view(set.zoom.limits = c(9,12))

There are many outer islands shown on the plot above. As seen from the plot, the areas to be excluded are: Pulau Harapan, Pulau Kelapa, Pulau Panggang, Pulau Tidung, Pulau Pari and Pulau Untung Jawa.

outer_islands <- c("PULAU HARAPAN", "PULAU KELAPA", "PULAU PANGGANG", "PULAU TIDUNG", "PULAU PARI", "PULAU UNTUNG JAWA")

Use the filter() function to remove the outer islands.

jkt <- jkt %>%
  filter(!DESA %in% outer_islands)

Now using the plot mode, plot out jkt and ensure all outer islands are excluded.

tmap_mode("plot")
tm_shape(jkt) +
  tm_polygons()

As you can see from the plot, all the outer islands are excluded.

2.2.3.2 Unnecessary Fields

Only the first nine columns of the sf dataframe is relevant. Using the select() function, pick out the first nine columns.

jkt <- select(jkt, 1:9)

2.2.4 Final jkt Data Frame

st_geometry(jkt)
Geometry set for 263 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -3644275 ymin: 663887.8 xmax: -3606237 ymax: 701380.1
Projected CRS: DGN95 / Indonesia TM-3 zone 54.1
First 5 geometries:

2.3 Riwayat File Vaksinasi DKI Jakarta

There are 12 excel files from July 2021 to June 2022 to be imported. Using read_excel() function, import the “Data Kelurahan” sheet.

july_2021 = read_excel("data/aspatial/july_2021.xlsx", sheet = 1)
aug_2021 = read_excel("data/aspatial/aug_2021.xlsx", sheet = "Data Kelurahan")
sept_2021 = read_excel("data/aspatial/sept_2021.xlsx", sheet = "Data Kelurahan")
oct_2021 = read_excel("data/aspatial/oct_2021.xlsx", sheet = "Data Kelurahan")
nov_2021 = read_excel("data/aspatial/nov_2021.xlsx", sheet = "Data Kelurahan")
dec_2021 = read_excel("data/aspatial/dec_2021.xlsx", sheet = "Data Kelurahan")
jan_2022 = read_excel("data/aspatial/jan_2022.xlsx", sheet = "Data Kelurahan")
feb_2022 = read_excel("data/aspatial/feb_2022.xlsx", sheet = "Data Kelurahan")
mar_2022 = read_excel("data/aspatial/mar_2022.xlsx", sheet = "Data Kelurahan")
apr_2022 = read_excel("data/aspatial/apr_2022.xlsx", sheet = "Data Kelurahan")
may_2022 = read_excel("data/aspatial/may_2022.xlsx", sheet = "Data Kelurahan")
june_2022 = read_excel("data/aspatial/june_2022.xlsx", sheet = "Data Kelurahan")

2.3.1 Data Cleaning

I will exclude all the outer_islands identified earlier.

july_2021 <- july_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
aug_2021 <- aug_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
sept_2021 <- sept_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
oct_2021 <- oct_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
nov_2021 <- nov_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
dec_2021 <- dec_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
jan_2022 <- jan_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
feb_2022 <- feb_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
mar_2022 <- mar_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
apr_2022 <- apr_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
may_2022 <- may_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
june_2022 <- june_2022 %>%
  filter(!KELURAHAN %in% outer_islands)

All the tables now have 262 observations, 1 less than the jkt multipolygon dataframe.

2.3.1.1 Unnecessary Columns

Only the first 6 columns are important to the exercise.

july_2021 <- select(july_2021, 1:6)
aug_2021 <- select(aug_2021, 1:6)
sept_2021 <- select(sept_2021, 1:6)
oct_2021 <- select(oct_2021, 1:6)
nov_2021 <- select(nov_2021, 1:6)
dec_2021 <- select(dec_2021, 1:6)
jan_2022 <- select(jan_2022, 1:6)
feb_2022 <- select(feb_2022, 1:6)
mar_2022 <- select(mar_2022, 1:6)
apr_2022 <- select(apr_2022, 1:6)
may_2022 <- select(may_2022, 1:6)
june_2022 <- select(june_2022, 1:6)

Now all the datasets only have 6 columns.

2.3.2 Combining Datasets

First, add a date column to each data frame using the mutate() and as.Date() function.

july_2021 <- july_2021 %>%
  mutate(dateMonth = as.Date("2021-07-01"))
aug_2021 <- aug_2021 %>%
  mutate(dateMonth = as.Date("2021-08-01"))
sept_2021 <- sept_2021 %>%
  mutate(dateMonth = as.Date("2021-09-01"))
oct_2021 <- oct_2021 %>%
  mutate(dateMonth = as.Date("2021-10-01"))
nov_2021 <- nov_2021 %>%
  mutate(dateMonth = as.Date("2021-11-01"))
dec_2021 <- dec_2021 %>%
  mutate(dateMonth = as.Date("2021-12-01"))
jan_2022 <- jan_2022 %>%
  mutate(dateMonth = as.Date("2022-01-01"))
feb_2022 <- feb_2022 %>%
  mutate(dateMonth = as.Date("2022-02-01"))
mar_2022 <- mar_2022 %>%
  mutate(dateMonth = as.Date("2022-03-02"))
apr_2022 <- apr_2022 %>%
  mutate(dateMonth = as.Date("2022-04-01"))
may_2022 <- may_2022 %>%
  mutate(dateMonth = as.Date("2022-05-01"))
june_2022 <- june_2022 %>%
  mutate(dateMonth = as.Date("2022-06-01"))

Next, I will combine all the 12 datasets into 1 dataset using rbind() function.

vacc_data <- rbind(july_2021, aug_2021, sept_2021, oct_2021, nov_2021, dec_2021, jan_2022, feb_2022, mar_2022, apr_2022, may_2022, june_2022)

2.3.3 Final Dataset

head(vacc_data)
# A tibble: 6 × 7
  `KODE KELURAHAN` `WILAYAH KOTA` KECAMATAN   KELUR…¹ SASARAN BELUM…² dateMonth 
  <chr>            <chr>          <chr>       <chr>     <dbl>   <dbl> <date>    
1 <NA>             <NA>           <NA>        TOTAL   7739060 5041111 2021-07-01
2 3172051003       JAKARTA UTARA  PADEMANGAN  ANCOL     20393   13272 2021-07-01
3 3173041007       JAKARTA BARAT  TAMBORA     ANGKE     25785   16477 2021-07-01
4 3175041005       JAKARTA TIMUR  KRAMAT JATI BALE K…   25158   18849 2021-07-01
5 3175031003       JAKARTA TIMUR  JATINEGARA  BALI M…    8683    5743 2021-07-01
6 3175101006       JAKARTA TIMUR  CIPAYUNG    BAMBU …   22768   15407 2021-07-01
# … with abbreviated variable names ¹​KELURAHAN, ²​`BELUM VAKSIN`

The final vacc_data dataset has a total of 7 columns:

  1. Identifiers: Kode Kelurahan / Kelurahan
  2. Targeted people to be vaccinated: Sasaran
  3. Yet to be vaccinated: Belum Vaksin
  4. Date identifier: dateMonth

3 Choropleth Mapping and Analysis

3.1 Monthly Vaccination Rates at Sub-district Level

3.1.1 Computing Monthly Vaccination Rates

  • SASARAN (Targeted people to be vaccinated)

  • BELUM VAKSIN (Yet to be vaccinated)

Using the equation below, compute the monthly vaccination rate and add the rate value as a new column ‘monthly_vacc_rate’ using the mutate() function.

\[ Monthly Vaccination Rate = \frac{(SASARAN - BELUMVAKSIN)}{SUM(SASARAN, BELUMVAKSIN)} \]

vacc_data <- vacc_data %>%
  mutate(monthly_vacc_rate = (SASARAN - `BELUM VAKSIN`)/(SASARAN + `BELUM VAKSIN`))

3.1.2 Combining Dataset and Shapefile

Using the left_join() function, with identifiers “DESA” and “KELURAHAN”, combine the vacc_data dataset with the jkt shapefile.

vacc_july_2021 <- vacc_data %>%
  filter(dateMonth == "2021-07-01")
jkt_july_2021 <- left_join(jkt, vacc_july_2021, by = c("DESA" = "KELURAHAN" ))

3.2 Choropleth Mapping

Using the qtm() function, plot out the choropleth map for July 2021.

qtm(jkt_july_2021, "monthly_vacc_rate", title = "July 2021", fill.palette="Greens")

As seen above, there are missing monthly vaccinated rate values for some sub-districts. There is a mismatch and lack of vaccination data for these areas.

I will leave them as missing values instead of giving a monthly vaccination rate of 0 to avoid confusion with actual areas with recorded 0 rates.

3.2.1 Other Months

Repeat the above steps for the next 11 months.

Code 
#AUGUST 2021
vacc_aug_2021 <- vacc_data %>%
  filter(dateMonth == "2021-08-01")
jkt_aug_2021 <- left_join(jkt, vacc_aug_2021, by = c("DESA" = "KELURAHAN" ))

#SEPTEMBER 2021
vacc_sept_2021 <- vacc_data %>%
  filter(dateMonth == "2021-09-01")
jkt_sept_2021 <- left_join(jkt, vacc_sept_2021, by = c("DESA" = "KELURAHAN" ))

#OCTOBER 2021
vacc_oct_2021 <- vacc_data %>%
  filter(dateMonth == "2021-10-01")
jkt_oct_2021 <- left_join(jkt, vacc_oct_2021, by = c("DESA" = "KELURAHAN" ))

#NOVEMBER 2021
vacc_nov_2021 <- vacc_data %>%
  filter(dateMonth == "2021-11-01")
jkt_nov_2021 <- left_join(jkt, vacc_nov_2021, by = c("DESA" = "KELURAHAN" ))

#DECEMBER 2021
vacc_dec_2021 <- vacc_data %>%
  filter(dateMonth == "2021-12-01")
jkt_dec_2021 <- left_join(jkt, vacc_dec_2021, by = c("DESA" = "KELURAHAN" ))

#JANUARY 2022
vacc_jan_2022 <- vacc_data %>%
  filter(dateMonth == "2022-01-01")
jkt_jan_2022 <- left_join(jkt, vacc_jan_2022, by = c("DESA" = "KELURAHAN" ))

#FEBRUARY 2022
vacc_feb_2022 <- vacc_data %>%
  filter(dateMonth == "2022-02-01")
jkt_feb_2022 <- left_join(jkt, vacc_feb_2022, by = c("DESA" = "KELURAHAN" ))

#MARCH 2022
vacc_mar_2022 <- vacc_data %>%
  filter(dateMonth == "2022-03-02")
jkt_mar_2022 <- left_join(jkt, vacc_mar_2022, by = c("DESA" = "KELURAHAN" ))

#APRIL 2022
vacc_apr_2022 <- vacc_data %>%
  filter(dateMonth == "2022-04-01")
jkt_apr_2022 <- left_join(jkt, vacc_apr_2022, by = c("DESA" = "KELURAHAN" ))

#MAY 2022
vacc_may_2022 <- vacc_data %>%
  filter(dateMonth == "2022-05-01")
jkt_may_2022 <- left_join(jkt, vacc_may_2022, by = c("DESA" = "KELURAHAN" ))

#JUNE 2022
vacc_june_2022 <- vacc_data %>%
  filter(dateMonth == "2022-06-01")
jkt_june_2022 <- left_join(jkt, vacc_june_2022, by = c("DESA" = "KELURAHAN" ))
Code 
qtm(jkt_july_2021, "monthly_vacc_rate", title = "July 2021", fill.palette="Greens")

Code 
qtm(jkt_aug_2021, "monthly_vacc_rate", title = "Aug 2021", fill.palette="Greens")

Code 
qtm(jkt_sept_2021, "monthly_vacc_rate", title = "Sept 2021", fill.palette="Greens")

Code 
qtm(jkt_oct_2021, "monthly_vacc_rate", title = "Oct 2021", fill.palette="Greens")

Code 
qtm(jkt_nov_2021, "monthly_vacc_rate", title = "Nov 2021", fill.palette="Greens")

Code 
qtm(jkt_dec_2021, "monthly_vacc_rate", title = "Dec 2021", fill.palette="Greens")

Code 
qtm(jkt_jan_2022, "monthly_vacc_rate", title = "Jan 2022", fill.palette="Greens")

Code 
qtm(jkt_feb_2022, "monthly_vacc_rate", title = "Feb 2022", fill.palette="Greens")

Code 
qtm(jkt_mar_2022, "monthly_vacc_rate", title = "Mar 2022", fill.palette="Greens")

Code 
qtm(jkt_apr_2022, "monthly_vacc_rate", title = "Apr 2022", fill.palette="Greens")

Code 
qtm(jkt_may_2022, "monthly_vacc_rate", title = "May 2022", fill.palette="Greens")

Code 
qtm(jkt_june_2022, "monthly_vacc_rate", title = "June 2022", fill.palette="Greens")

3.3 Spatial Patterns

Over the year, the monthly vaccination rate has been on an overall increase as seen in the legend where the minimum rate rises from 0.10 in July 2021 to 0.6 in June 2022.

Sub-district with consistently high monthly vaccination rates (relative to other areas) is Kelapa Gading Timur. Whereas the sub-district with consistently low monthly rates (relative to other areas) includes Kebon Melati and Petamburan.

March 2022 had the highest monthly vaccination rates.

Overall, excluding the areas with missing values, by October 2021, the monthly vaccination rate across Jakarta has been meeting half of the targeted numbers and has only been increasing since.

4 Local Gi* Analysis

4.1 Computing Local Gi* of Monthly Vaccination Rate

To compute local Gi* statistics, it requires the neighbour list (using st_contiguity()) and the wts (using st_inverse_distance()).

Using jkt_july_2021 as an example:

jkt_july_2021 = jkt_july_2021 %>% drop_na(monthly_vacc_rate)
jkt_july_2021_idw <- jkt_july_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )
Note

There are ‘NA’ values in the monthly_vacc_rate column due to missing Kelurahans of the two datasets. These rows need to be dropped to compute the nb and wts.

Next, calculate the local Gi* values using local_gstar_perm() with the nb and wt values computed in the previous step. Add the local Gi* values into a new sf dataset HCSA_july_2021 as a new column with the mutate() function, adding .before = 1 to place the new column as the first column. The unnest() function is used to make each element of the local_Gi list its own row.

HCSA_july_2021 <- jkt_july_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)
HCSA_july_2021
Simple feature collection with 252 features and 26 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -3644275 ymin: 663887.8 xmax: -3606237 ymax: 701380.1
Projected CRS: DGN95 / Indonesia TM-3 zone 54.1
# A tibble: 252 × 27
   gi_star    e_gi      var_gi p_value p_sim p_fol…¹ skewn…² kurto…³ nb    wts  
     <dbl>   <dbl>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>   <dbl> <nb>  <lis>
 1   2.08  0.00393     1.04e-7 3.73e-2  0.08    0.04  0.0725  -0.161 <int> <dbl>
 2   3.99  0.00395     5.02e-8 6.47e-5  0.02    0.01  0.371   -0.563 <int> <dbl>
 3  -0.496 0.00396     9.47e-8 6.20e-1  0.68    0.34  0.164   -0.463 <int> <dbl>
 4  -0.906 0.00396     1.09e-7 3.65e-1  0.44    0.22  0.138   -0.826 <int> <dbl>
 5   2.51  0.00396     5.67e-8 1.22e-2  0.02    0.01  0.192   -0.584 <int> <dbl>
 6   1.52  0.00395     6.58e-8 1.29e-1  0.16    0.08  0.0705  -0.216 <int> <dbl>
 7   0.963 0.00396     7.06e-8 3.36e-1  0.36    0.18 -0.106   -0.572 <int> <dbl>
 8  -0.125 0.00400     1.11e-7 9.00e-1  0.9     0.45  0.0655  -0.169 <int> <dbl>
 9   0.714 0.00400     6.90e-8 4.75e-1  0.5     0.25  0.354   -0.395 <int> <dbl>
10  -0.306 0.00397     8.26e-8 7.60e-1  0.7     0.35 -0.0967  -0.521 <int> <dbl>
# … with 242 more rows, 17 more variables: OBJECT_ID <dbl>, KODE_DESA <chr>,
#   DESA <chr>, KODE <dbl>, PROVINSI <chr>, KAB_KOTA <chr>, KECAMATAN.x <chr>,
#   DESA_KELUR <chr>, JUMLAH_PEN <dbl>, `KODE KELURAHAN` <chr>,
#   `WILAYAH KOTA` <chr>, KECAMATAN.y <chr>, SASARAN <dbl>,
#   `BELUM VAKSIN` <dbl>, dateMonth <date>, monthly_vacc_rate <dbl>,
#   geometry <MULTIPOLYGON [m]>, and abbreviated variable names ¹​p_folded_sim,
#   ²​skewness, ³​kurtosis

4.2 Displaying Gi* Maps

Using the plot mode, plot out both the Gi* Map and the p-value of Gi* using the HCDA_july_2021 sf dataset containing the newly added Gi* values. I will be using the “-RdBu” palette to better visualise the extreme positive and negative Gi* values which represents the hot and cold spot clusters respectively.

tmap_mode("plot")
map1 <- tm_shape(HCSA_july_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of July 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_july_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3 Other Months

Repeat the previous steps to the other 11 months.

4.3.1 August 2021

Code 
jkt_aug_2021 = jkt_aug_2021 %>% drop_na(monthly_vacc_rate)
jkt_aug_2021_idw <- jkt_aug_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_aug_2021 <- jkt_aug_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_aug_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of August 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_aug_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.2 September 2021

Code 
jkt_sept_2021 = jkt_sept_2021 %>% drop_na(monthly_vacc_rate)
jkt_sept_2021_idw <- jkt_sept_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_sept_2021 <- jkt_sept_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_sept_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of September 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_sept_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.3 October 2021

Code 
jkt_oct_2021 = jkt_oct_2021 %>% drop_na(monthly_vacc_rate)
jkt_oct_2021_idw <- jkt_oct_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_oct_2021 <- jkt_oct_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_oct_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of October 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_oct_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.4 November 2021

Code 
jkt_nov_2021 = jkt_nov_2021 %>% drop_na(monthly_vacc_rate)
jkt_nov_2021_idw <- jkt_nov_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_nov_2021 <- jkt_nov_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_nov_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of November 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_nov_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.5 December 2021

Code 
jkt_dec_2021 = jkt_dec_2021 %>% drop_na(monthly_vacc_rate)
jkt_dec_2021_idw <- jkt_dec_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_dec_2021 <- jkt_dec_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_dec_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of December 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_dec_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.6 January 2022

Code 
jkt_jan_2022 = jkt_jan_2022 %>% drop_na(monthly_vacc_rate)
jkt_jan_2022_idw <- jkt_jan_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_jan_2022 <- jkt_jan_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_jan_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of January 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_jan_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.7 February 2022

Code 
jkt_feb_2022 = jkt_feb_2022 %>% drop_na(monthly_vacc_rate)
jkt_feb_2022_idw <- jkt_feb_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_feb_2022 <- jkt_feb_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_feb_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of February 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_feb_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.8 March 2022

Code 
jkt_mar_2022 = jkt_mar_2022 %>% drop_na(monthly_vacc_rate)
jkt_mar_2022_idw <- jkt_mar_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_mar_2022 <- jkt_mar_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_mar_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of March 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_mar_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.9 April 2022

Code 
jkt_apr_2022 = jkt_apr_2022 %>% drop_na(monthly_vacc_rate)
jkt_apr_2022_idw <- jkt_apr_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_apr_2022 <- jkt_apr_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_apr_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of April 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_apr_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.10 May 2022

Code 
jkt_may_2022 = jkt_may_2022 %>% drop_na(monthly_vacc_rate)
jkt_may_2022_idw <- jkt_may_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_may_2022 <- jkt_may_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_may_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of May 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_may_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.11 June 2022

Code 
jkt_june_2022 = jkt_june_2022 %>% drop_na(monthly_vacc_rate)
jkt_june_2022_idw <- jkt_june_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_june_2022 <- jkt_june_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_june_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "Gi* of June 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_june_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.4 Statistical Conclusion

At the start of the year, the hot and cold spots were spreaded across the sub-districts. But overtime, the Gi* value bottom half of Jakarta has been significantly more positive. More specifically areas such as Jagakarsa, Ciganjur and their neighbours have been a extremely significant positive cluster especially at the later half of the year. This means that these sub-districts contains the highest monthly vaccination rates as a cluster as compared to other sub-district areas.

Kebon Melati has been a significant cold spot consistently throughout the year. Its neighbouring sub-districts are also visualised as cold spots throughout the year but is less significant than Kebon Melati itself. This means that Kebon Melati and its surrounding neighbours contains significantly lower monthly vaccination rates as a cluster as compared to other sub-district areas.

5 Emerging Hot Spot Analysis

5.1 Mann-Kendall Test

To perform the Mann-Kendall test, a time series cube needs to be created first.

5.1.1 Creating Time Series Cube

Since the spacetime() function only accepts the location column of 2 datasets to be the same, I will update the jkt “KESA” column to “KELURAHAN” for standardisation using the colnames() function.

colnames(jkt)[3] <- "KELURAHAN"

Next, using the spacetime() function, create a spacetime object with “KELURAHAN” as the location column and “dateMonth” as the time column.

#vacc_data = vacc_data %>% drop_na(c("KODE KELURAHAN", "WILAYAH KOTA", "KECAMATAN"))
#vacc_data_st <- spacetime(vacc_data, jkt, ".KELURAHAN", "dateMonth")
vacc_data_st <- spacetime(
  .data = vacc_data, 
  .geometry = jkt, 
  .loc_col = "KELURAHAN", 
  .time_col = "dateMonth")

There are missing values in the spacetime object. It will not be a spacetime cube due to the lack of rows. Using complete_spacetime_cube() function, create a completed spacetime cube vacc_data_st_complete.

vacc_data_st_complete <- complete_spacetime_cube(vacc_data_st)
is_spacetime_cube(vacc_data_st_complete)
[1] TRUE

Now the spacetime cube is valid.

5.1.2 Computing Spatial Weights

Before implementing the Mann-Kendall test, we need the nb and wt which will be computed similarly to the method use in the previous sections.

vacc_st_nb <- vacc_data_st_complete %>%
  activate("geometry") %>%
  mutate(nb = include_self(st_contiguity(geometry)),
         wt = st_inverse_distance(nb, geometry,
                                  scale = 1,
                                  alpha = 1),
         .before = 1) %>%
  set_nbs("nb") %>%
  set_wts("wt")
head(vacc_st_nb)
# A tibble: 6 × 10
  KELUR…¹ dateMonth  KODE …² WILAY…³ KECAM…⁴ SASARAN BELUM…⁵ month…⁶ nb    wt   
  <chr>   <date>     <chr>   <chr>   <chr>     <dbl>   <dbl>   <dbl> <lis> <lis>
1 KEAGUN… 2021-07-01 317303… JAKART… TAMAN …   15262   10275   0.195 <int> <dbl>
2 GLODOK  2021-07-01 317303… JAKART… TAMAN …    7192    3561   0.338 <int> <dbl>
3 HARAPA… 2021-07-01 317103… JAKART… KEMAYO…   19987   13623   0.189 <int> <dbl>
4 CEMPAK… 2021-07-01 317103… JAKART… KEMAYO…   27330   18417   0.195 <int> <dbl>
5 PASAR … 2021-07-01 317102… JAKART… SAWAH …   11656    6067   0.315 <int> <dbl>
6 KARANG… 2021-07-01 317102… JAKART… SAWAH …   23056   15040   0.210 <int> <dbl>
# … with abbreviated variable names ¹​KELURAHAN, ²​`KODE KELURAHAN`,
#   ³​`WILAYAH KOTA`, ⁴​KECAMATAN, ⁵​`BELUM VAKSIN`, ⁶​monthly_vacc_rate

5.1.3 Computing Gi*

Next, I will calculate the Gi* value and add it into a new dataframe called gi_stars. Since there are missing values in the monthly_vacc_rate, I will drop the empty rows before calulating the Gi* values with the same method as the previous section.

vacc_st_nb <- vacc_st_nb %>% drop_na("monthly_vacc_rate")
gi_stars <- vacc_st_nb %>% 
  mutate(gi_star = local_gstar_perm(
    monthly_vacc_rate, nb, wt)) %>% 
  unnest(gi_star)

5.1.4 Mann-Kendall Test

To perform the Mann-Kendall test, I will call the MannKendall() function on the newly calculated Gi* values into a list mk. I will then add this computed results into a new dataframe ehsa.

ehsa <- gi_stars %>%
  group_by(KELURAHAN) %>%
  summarise(mk = list(
    unclass(MannKendall(gi_star)))) %>%
  tidyr::unnest_wider(mk)

5.3 EHSA Map

Lastly, I will run a emerging hotspot analysis using emerging_hotspot_analysis() function on the “monthly_vacc_rate” variable of the completed spacetime cube, running for 100 tests.

vacc_data_st_complete$monthly_vacc_rate[is.na(vacc_data_st_complete$monthly_vacc_rate)] <- 0
ehsa <- emerging_hotspot_analysis(
  x = vacc_data_st_complete, 
  .var = "monthly_vacc_rate", 
  k = 1, 
  nsim = 99
)

I will then combine the computed ehsa into the jkt shapefile using the left_join() function with columns “KELURAHAN” from jkt and “location” from ehsa.

jkt_ehsa <- jkt %>%
  left_join(ehsa, by = c("KELURAHAN" = "location" ))

To plot out the significant emerging hotspot analysis results, I will filter out the results that have p-values lesser than 0.05 with the filter() function.

ehsa_sig <- jkt_ehsa  %>%
  filter(p_value < 0.05)
tmap_mode("view")
tm_shape(jkt_ehsa) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
  tm_text("KELURAHAN", size = 0.7) +
tm_shape(ehsa_sig) +
  tm_fill("classification") + 
  tm_borders(alpha = 0.4) +
  tm_view(set.zoom.limits = c(11,14))

5.4 Spatial Patterns

Majority of the sub-districts have a classification of oscillating hotspot. This aligns with the previous Local Gi* maps with the bottom half being majority of the hotspot.

Some notable sub-districts include:

  • Petojo Utara: Persistent coldspot

  • Pademangan Timur: Intensifying hotspot

It is interesting how the areas I have mentioned to be significant earlier not remotely near the 2 sub-districts computed to be significant here.

6 Conclusion

For each test, there were different sub-districts that I perceived or was computed to be significant. It was different for every single test which made me realise the importance of doing multiple tests to have a better understanding of the data I am handling.